# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from hysop.tools.htypes import first_not_None, check_instance
from hysop.tools.sympy_utils import greak, Greak, subscripts
from hysop.fields.continuous_field import Field, TensorField
from hysop.tools.numpywrappers import npw
from hysop.constants import BoxBoundaryCondition, BoundaryCondition
[docs]
def VelocityField(domain, name=None, pretty_name=None, idx=None, **kwds):
name = first_not_None(name, "U")
pretty_name = first_not_None(pretty_name, greak[20])
name += "" if (idx is None) else str(idx)
pretty_name += "" if (idx is None) else f"[{idx}]"
lboundaries, rboundaries = domain.lboundaries, domain.rboundaries
dim = domain.dim
def velocity_boundaries(boundaries, component):
check_instance(boundaries, npw.ndarray, values=BoxBoundaryCondition)
fboundaries = npw.empty_like(boundaries)
fboundaries[...] = None
for i, bd in enumerate(boundaries):
is_normal = dim - i - 1 == component
if bd is BoxBoundaryCondition.PERIODIC:
fbd = BoundaryCondition.PERIODIC
elif (
bd is BoxBoundaryCondition.SYMMETRIC
): # (normal to boundary velocity = 0)
if is_normal:
fbd = BoundaryCondition.HOMOGENEOUS_DIRICHLET
else:
fbd = BoundaryCondition.HOMOGENEOUS_NEUMANN
elif bd is BoxBoundaryCondition.OUTFLOW: # (velocity is normal to boundary)
if is_normal:
fbd = BoundaryCondition.HOMOGENEOUS_NEUMANN
else:
fbd = BoundaryCondition.HOMOGENEOUS_DIRICHLET
else:
msg = "FATAL ERROR: Unknown domain boundary condition {}."
msg = msg.format(bd)
raise NotImplementedError(msg)
fboundaries[i] = fbd
check_instance(fboundaries, npw.ndarray, values=BoundaryCondition)
return fboundaries
def _make_field(idx, **fkwds):
# Adapt velocity boundaries to domain boundaries
(component,) = idx
fkwds["lboundaries"] = velocity_boundaries(lboundaries, component)
fkwds["rboundaries"] = velocity_boundaries(rboundaries, component)
return TensorField.default_make_field(idx=idx, **fkwds)
kwds.setdefault("make_field", _make_field)
kwds.setdefault("is_vector", True)
return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
[docs]
def VorticityField(velocity, name=None, pretty_name=None, idx=None, **kwds):
# vorticity domain domain boundaries are deduced from velocity boundary conditions
check_instance(velocity, Field)
domain = velocity.domain
assert velocity.nb_components == domain.dim, "Invalid velocity Field."
name = first_not_None(name, "W")
pretty_name = first_not_None(pretty_name, greak[24])
name += "" if (idx is None) else str(idx)
pretty_name += "" if (idx is None) else f"[{idx}]"
return velocity.curl(
name=name,
pretty_name=pretty_name,
scalar_name_prefix=name,
scalar_pretty_name_prefix=pretty_name,
**kwds,
)
[docs]
def DensityField(domain, name=None, pretty_name=None, **kwds):
name = first_not_None(name, "rho")
pretty_name = first_not_None(pretty_name, greak[16])
return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
[docs]
def ViscosityField(domain, name=None, pretty_name=None, mu=False, **kwds):
if mu:
name = first_not_None(name, "mu")
pretty_name = first_not_None(pretty_name, greak[11])
else:
name = first_not_None(name, "nu")
pretty_name = first_not_None(pretty_name, greak[12])
return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
[docs]
def LevelSetField(domain, name=None, pretty_name=None, **kwds):
name = first_not_None(name, "phi")
pretty_name = first_not_None(pretty_name, Greak[21])
return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
[docs]
def PenalizationField(domain, name=None, pretty_name=None, **kwds):
name = first_not_None(name, "lambda")
pretty_name = first_not_None(pretty_name, greak[10])
return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
[docs]
def CurvatureField(domain, name=None, pretty_name=None, **kwds):
name = first_not_None(name, "kappa")
pretty_name = first_not_None(pretty_name, greak[9])
return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)